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We calculate the first six non-normalized moments of particle multiplicity within the hadron 
resonance gas model. In terms of lower moments, general expressions for arbitrary higher moments 
are deduced, so that it is concluded that going from lower to higher moments is achievable through a 
series of lower moments and correlation functions. The thermal evolution of the first four normalized 
moments and their products (ratios) are studied at different chemical potentials /i, so that it was 
possible to evaluate them at the chemical freeze-out curve, which is characterized by constant s/T 3 
■ at all values of n, where s and T are entropy density and temperature, respectively. It has been 

found that non-monotonic behavior reflecting dynamical fluctuation and strong correlations appears 
starting from the normalized third moment (skewness S). The signatures of non-monotonicity in 
00 ' the normalized fourth moment (kurtosis re) and its products are very clear. Based on these findings, 

we introduce novel conditions describing the chemical freeze-out curve and the location of the QCD 
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critical endpoint. The chemical freeze-out curve is described by re = 3 \ 2 > where \ is the susceptibility 
in particle number. The location of QCD critical endpoint (at T- and /x-axis) is found when Sb = Sf 
or Kb = re/, where the subscripts b and / refer to bosons and fermions, respectively. Accordingly, 
Q y the QCD endpoint is positioned at fj, ~ 350 MeV and T ~ 162 MeV. 

r~| ■ PACS numbers: 74.40.Gh, 05.70.Fh 

Keywords: Fluctuation phenomena non-equilibrium processes, Phase transitions in statistical mechanics and 
thermodynamics 



l>- . I. INTRODUCTION 

The meanings of "moments" are rooted in essential mathematical aspects. Physics provides with accurate 
expressions for the "moments". Based on this, we find that "moments" have different definitions in physics 
and mathematics. In mathematics, they are the quantitative estimation for the potential shape of a certain set 
of measurements. For example, the widely used "first and second moments" measure the "distribution" and 
the "width" of a set of measurements in various dimensions, respectively. Obviously, the higher moments are 
supposed to cover other aspects. For instance, the way that the distribution is skewed from its mean, or peaked 
around it, can be described by "third" and "fourth moments", respectively. In physics, the "moment" is often 
represented in a somewhat different way. In general, any distribution can be characterized by the moments of a 
function describing the nature of the distribution itself and a number of other features, such as mean, variance, 
skewness, etc. It is obvious that the higher the moments (especially the ones beyond the fourth moment), the 
harder to estimate them. It is always valid that large samples are required in order to estimate them, correctly. 

Recently, the higher moments have gained prominence in high energy physics with a huge hope in pinpointing 
the QCD critical endpoint (CEP) [l[ connecting the first-order boundary separating the hadronic from the 
partonic matter at high density with the cross-over boundary at low density 0, H- There are large theoretical 
uncertainties of its location and even its existence is not fully confirmed, yet [4|-|6| . The characteristic feature 
of CEP are critical dynamical fluctuations @, The higher moments are conjectured to reflect the large 
fluctuations associating the hadron-quark phase transition. This was the motivation of a remarkable number 
of experimental and theoretical studies [l|, l9rlll|. Recently, various calculations have shown that the higher 
moments of the multiplicity distributions of some conserved quantities, such as net-baryon, net-charge, and 
net-strangeness, are sensitive to the correlation length £ [l2l - fl4j . which in turn is related to the higher moments, 
themselves. In realistic heavy-ion collisions, the correlation length is found to remain finite. 
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Apparently, it should not be an exception that the strongly interacting QCD matter undergoes phase transi- 
tion^) at extreme conditions, as almost all matter types suffer from such a critical change as the extreme con- 
ditions change. The Lattice QCD calculations predict that a cross-over takes place between the hadronic phase 
and the Quark Gluon Plasma (QGP) phase, when the temperature exceeds critical value of T c ~ 150—190 MeV. 
As given in Ref. HEi, lattice Q CD 

assigns different values to T c depending on different parameters (for in- 
stance, the quark favours and their masses) and partly the lattice configurations. Furthermore, the values of 
T c depend on the baryon chemical potential (i. At vanishing chemical potential, the lattice QCD calculations 
show that the higher order susceptibilities can be related to the higher order moments of the corresponding mul- 
tiplicity distributions. Furthermore, they show non-monotonic behavior near the critical temperature T c 
Apparently, the QCD phase diagram including CEP can be mapped out using this characteristic behaviour. In 
literature, one finds that the QCD phase diagram has been determined using various criteria, such as critical 
energy density [lTj-Qjl or line of constant physics Q . Also the chemical freeze-out curve can be characterized 
assuming constant physics. For a recent and extensive review, we refer to [2(| and the references therein. 

At large chemical potential fi, various calculations based on QCD-based models indicate that the transition 
from the hadronic phase to the QGP phase is of first-order. The endpoint connecting this line with the one 
of the cross over (at small /x) is likely of second-order 0, H, HH H3|- Therefore, it should be accompanied with 
large fluctuations in different physical quantities. Experimentally, we study the QCD phase diagram by varying 
the colliding energy in heavy- ion collisions and determining the critical temperature T c [24| . The possibility of 
finding QCD CEP has motivated the RHIC beam energy scan program [25[. By tuning the collision energy from 
a center-of-mass energy of 200 GeV down to 5 GeV, we are able to vary the baryon chemical potential from 
~ 2 to ~ 500 MeV. Nevertheless, there are large theoretical uncertainties of its location and even its existence is 
not fully confirmed [1, H, Hi], [23| . The higher moments of various physical quantities are conjectured to reflect 
the non-monotonic behaviour near the QCD CEP. In present work, we give a systematic study for the higher 
moments of the multiplicity distribution and show their ability in reflecting the fluctuations along the QCD 
phase diagram. 

In this paper, we study the first six non-normalized higher moments of the particle multiplicity in the hadron 
resonance gas (HRG) model, which is introduce in section [TTJ In section III Ai we suggest an answer to the 
question: What is needed when going from "trivial" lower to "sophisticated" higher moments? Based on 
results, general expressions for arbitrary higher moments are deduced, so that it is possible to conclude that 
going from lower to higher moments is achievable through a series of lower moments and correlation functions. 
The thermal evolution of the first four normalized moments and their products (ratios) are studied at different 
chemical potentials fi in section Hi Bl Therefore, it is possible to evaluate them at the chemical freeze-out curve, 
which is characterized by constant s/T 3 at all values of //, where s and T are entropy density and temperature, 
respectively Sections IIIII and IIVI are devoted to introduce novel conditions describing the chemical freeze-out 
curve and define the location of the QCD critical endpoint. 



II. HIGHER MOMENTS OF THE PARTICLE MULTIPLICITY 



In a grand canonical ensemble, it is straightforward to derive an expression for the free energy. In section 
III Al we give a list of different moments of particle multiplicity. It seems that going from trivial "lower" 
to sophisticated "higher" orders is accompanied by adding/subtracting certain distribution functions. These 
functions seem to reflect the higher order correlations. In section IIIBI we list out various normalized higher 
moments, where the normalization is based on the variance a. The products (ratios) of normalized moments 
are elaborated in section Hi CI 

The hadron resonances treated as a free gas [H, [l7l - [l9l [2r| are conjectured to add to the thermodynamic 
pressure in the hadronic phase (below T c ). This statement is valid for free as well as strong interactions between 
the resonances themselves. It has been shown that the thermodynamics of strongly interacting system can also 
be approximated to an ideal gas composed of hadron resonances with masses < 2 GeV @, H3|. Therefore, the 
confined phase of QCD, the hadronic phase, is modelled as a non-interacting gas of resonances. The grand 
canonical partition function reads 

Z(T, V) = Tr [cxp- H/T ] (1) 

where H is the Hamiltonian of the system and T is the temperature. The Hamiltonian is given by the sum of 
the kinetic energies of relativistic Fermi and Bose particles. The main motivation of using this Hamiltonian is 
that it contains all relevant degrees of freedom of confined and strongly interacting matter. It includes implicitly 
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the interactions that result in resonance formation. In addition, it has been shown that this model can submit 
a quite satisfactory description of particle production in heavy-ion collisions. With the above assumptions 
the dynamics the partition function can be calculated exactly and be expressed as a sum over single-particle 
partition functions Z\ of all hadrons and their resonances. 

\nZ(T,fa,V) = J2^Zl(T,V)=J2±^ f°° k 2 dkln{l±exp[( fH -e i )/T}} (2) 

i i 

where ej(fc) = (k 2 + to?) 1 / 2 is the i— th particle dispersion relation, gi is spin-isospin degeneracy factor and ± 
stands for bosons and fcrmions, respectively. 

Before the discovery of QCD, it was speculated about a possiblephase transition of a massless pion gas to a 
new phase of matter. Based on statistical models like Hagedorn [28| and Bootstrap [2{|, the thermodynamics 
of such an ideal pion gas is studied, extensively. After the QCD, the new phase of matter is known as quark 
gluon plasma (QGP). The physical picture was that at T c the additional degrees of freedom carried by QGP 
are to be released resulting in an increase in the thermodynamic quantities like energy and pressure densities. 
The success of hadron resonance gas model in reproducing lattice QCD results at various quark flavours and 
masses (below T c ) changed this physical picture drastically. Instead of releasing additional degrees of freedom 
at T > T c , the interacting system reduces its effective degrees of freedom at T < T c . In other word, the hadron 
gas has much more degrees of freedom than QGP. 

At finite temperature T and baryo-chemical potential fii , the pressure of i-th hadron or resonance reads 

poo 

p(T,m) = ±f,T k 2 dk\n{l±exp[(fa-e l )/T}}. (3) 
27r Jo 

As no phase transition is conjectured in HRG, summing over all hadron resonances results in the final thermo- 
dynamic pressure in the hadronic phase. 

The switching between hadron and quark chemistry is given by the relations between the hadronic chemical 
potentials and the quark constituents; fa = Znt>pL q + n s [j,g, where rib{n s ) being baryon (strange) quantum 
number. The chemical potential assigned to the light quarks is fi q = (/j, u + Hd)/2 and the one assigned to 
strange quark reads ^5 = fj, q — fx s . The strangeness chemical potential ^5 is calculated as a function of T 
and fa under the assumption that the overall strange quantum number has to remain conserved in heavy-ion 
collisions @. 



A. Non-Normalized higher moments 

As given in the introduction, the higher moments can be studied in different physical quantities. For example, 
the higher moments of charged-particle multiplicity distribution have been predicted four decades ago [30j . 
Recently, the higher moments of various multiplicity distributions have been reported by STAR measurements 
[9l. rToj| and lattice QCD calculations (3l|- [33 . The empirical relevance of the higher moments to the experimental 
measurements has been suggested in Ref . [l| . Accordingly, the measurement of the correlation length seems to be 
very much crucial. In a future work, we shall discuss the experimental signatures of the higher moments. Another 
ingredient would be the experimental sensitivity for the suggested signatures, which are based on singular 
behavior of thermodynamic functions. For the time being, we just want to mention that the experimental 
measurements apparently take place at the final state of the collision, which would mean that the signals have 
to survive the extreme conditions in such high energy collisions. It has been pointed out that the contribution of 
the critical fluctuations to the higher moments is proportional to a positive power of £. The latter is conjectured 
to diverge at QCD CEP and such an assumption is valid in the thermodynamic limit. 

In present work, we first study the contribution of the hadron resonance mass spectrum in clarifying the 
interrelations between the various non-normalized and normalized moments. The motivation of using resonance 
gas is clear fl7rfl9l ]. As introduced about five decades ago, the resonance gas is only valid below T c 0, Il7l - 
GiiHE HI]. Therefore, present work is valid in the hadronic phase, only. 

For z-th particle, the " first" moment is given by the derivative of p = —Td In Z(T, V,fa)/ dV with respect to 
the dimensionless quantity fa/T. When taking into account the antiparticles, we add a negative sign to the 
chemical potential. The first derivative describes the multiplicity distribution or an expectation operator, which 
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is utilized to estimate the number or multiplicity density 

ft rr r ^ k2dk 



f c 

mi(T,Hi) = ±7r^ T / — 
2 ?r 2 Jo 1 



± e t 

The "second" moment is known as the variance. It gives the susceptibility of the measurements. 

g l m f°° k 2 dk _ g, „ r e 2 ^ k 2 dk 

lo \±e ! ^ 2lr " Jo fiie^ 



(4) 



The "third" moment measures of the lopsidedness of the distribution. As given in section lllBl the normalization 
of third moment is known as skewness or the asymmetry in the distribution. Skewness tells us about the direction 
of variation of the data set. 



x , ft ™ f°° e ^ edk 9i o T f°° e 2 ^ k 2 dk 9i [°° e 3 ^ k 2 dk 



l±e^ ™ •>* l±eV 



In general, the normalization of r— th moment is obtained by dividing it by o~ r , where cr is the standard deviation. 
The normalization is assumed to remove the brightness dependence. 

The "fourth" moment compares the tallness and skinny or shortness and squatness, i.e. shape, of a certain 
measurement to its normal distribution. It is defined as the uncertainty is an uncertainty or "the location- and 
scale-free movement of probability mass from the shoulders of a distribution into its center and tails and to 
recognize that it can be formalized in many ways " [3^ . 

9i /'°° k 2 dk 9l [°° e 2 ^ k 2 dk g [°° e 3 ^ k 2 dk 



o (l±e^) 2?r ^ (l±e^) 



2 6T / —j. (7) 



ft 

l±e t 

The "fifth" moment measures the asymmetry sensitivity of the "fourth" moment. 

e ^ k dk gi I e t k dk . g% e t k dk 



trp \ j. 9i rr f e T k dk ft KT /* e t k 2 dk g, f 
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The "sixth" moment is generally associated with compound options. 

g p e^fc 2 rffc <?; /'°° e 2 ^ k 2 dk 9l [°° e 3 ^ PA 



ft ,„ OT A 00 e4 ^" fc2dfc ^ ft /'°° e 5 ^" k 2 dk 9l 1<MVf , f°° e 6 ^ fc 2 dfc 

1 ± e -1 ' 7 



9 -390T / j ± -^360T / % - -^120T / (9) 

2- 2 io (l±e^) 4 ^ ^ (l±e^) 5 2 ^ 



Thus, from Eqs. ((U)-©, a general expression for r-th moment can be deduced 

i T Y, arl / T 
i=i Jo ( i 



m r (T,fXi) = ^a r .i j — t k z dk, (10) 



± e T 
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where the coefficients read 



a r ,i = (±l)'(-l) ,+1 pa r _i 1 i + (l-l)o P _i,i_ 1 ] 



(11) 



where I < r and a r j vanishes Vr < 1. It is obvious that the coefficients of a certain moment are to be determined 
from a long chain of all previous ones. Such a conclusion would date back to about four decades [3(|. It has been 
shown that the coefficients are related to high order correlation functions. Should this assumption is proven to 
be valid, then expression (1111) gets a novel interpretation. It seems to sum up the correlation functions up to the 
7'-th order. According to [30( and when neglecting the two-particle correlations C2, then the higher moments 
read 



m 2 
m 3 

J7l4 



((6N) 2 ) 
((SNf) 
((5Nf) 



«4(A)+C 3 , 

= 6m 3 + 3m 2 . — 8m-2 + 8C4 



where SN = A - ( N) , 



C 3 (Pl,P2,P3) = ^2(PlPlPl)+3 X] (PlP2P2)+6 2J (P1P2P3), 

Pi Pl<P2 Pl<P2<P3 

Ci(pi,P2,P3,P4,) = ^2 (P1P2P2P2) + + 8 ^ (PlP2P3P4), 



(12) 
(13) 
(14) 



(15) 
(16) 



Pl<P2 



Pl<P2<P3 



Pl<P2<P3<P4 



and pi is the i—th particle. To the second moment m-i we have to add the effects of the two particle correlation 
function, 2 Ci- The third moment m 3 gets approximately three times this amount. The three and four particle 
correlations, Eqs. (|15[) and (|16|) . appear first in third and fourth moment, Eq. (|14[) and (|15p . respectively. 

When ignoring these high order correlations, we are left with the widely used two particle correlations. In 
present work, we restrict the discussion to this type, only. It is of great interest, as it is accessible experimentally 
and achievable, numerically. The two particle correlations are suggested as a probe for the bulk QCD medium, 
energy loss, medium response, jet properties and intensity interferometry [35j-|4l|. In addition to this list of 
literature, the comprehensive review |42J can be recommended. 
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Fig. 1: Left panel depicts dimensionless non-normalized moments in dependence on the temperature. Right panel shows 
the evolution of the so-called correlation parts (the integrals appearing in Eqs. (|20[) -(125 [ )) with the temperature. 



Taking into account the particle multiplicities, then expressions (JT 
follows. 



77J4 



m 3 - ((MO 2 ) 
m 3 = ((SNf) 
((5N) 4 ) - 3((6N) 2 } 2 



(N 2 ) - (AO 2 , 

(A 3 ) - (A 2 ) (A) +2(A 3 ) 3 , 

((A) 4 ) - 2((A 2 ) 2 ) - 5((A 2 )) 2 + 6((A) 2 ) ((A 2 )). 



T5|) and (fbfj) can be re-written as 



(17) 
(18) 
(19) 



As mentioned in section III Al the experimental signals for the higher moments shall be discussed in a future 
work. Another conclusion can be made so far. It is obvious that expressions (|20)) - (f2"5"|) lead to the conclusion 



6 



that all moments are entirely depending on the previous ones. Thus, the r-th moment has contributions coming 
from all moments with orders lower than r. 

As given in Eq. (jlOp . the dependence of first six non- normalized moments on lower ones can explicitly be 
deduced as follows. 



fli 



e ? k dk 



mr = j^T I — , (20) 



o l±e ? 
9i 



f°° e 2 ^ k 2 dk 



lie t 

m 3 = -2m, + 3m 2 + 2^T / - f (22) 

27H ./ , fs^iiV 



lie t 



6mi — llm 2 + 6m 3 — 6— -^-T / — 



7 7T' 

± e l ~ 



m 5 = -24mi + 50m 2 - 35m 3 + 10m 4 + / ^-^f, (24) 

2 ^ 7 ° ^lie^ 



oo o5 ^±ii fc 2 dfc 



. r°° e 6 jfc 2 dA; 
m 6 = 120mi - 274m 2 + 225m 3 - 85m 4 + 15mr 3 - 120^rT / (25) 

27T J ° fi± e :£ ^y 



Naively spoken, we conclude that raising lower moments to higher ones is achievable through a series of all 
lower moments and an additional term reflecting the correlations, themselves [30]. From Eqs. (P0)) - ([23|) . the 
additional terms are proportional to (N) r . Apparently, they are essential in order to judge whether going from 
lower to higher moments would make the signatures of dynamical fluctuations clearer than when excluding 
them. 

The first and second terms can be generalized. Then, a general expression would read 

m r = (-) r " 1 (r - ly-^mx - [(r - l)cC + (r - 2)!] m 2 + 0(m >2 ), (26) 

where c mr is the coefficient of r-th moment. The last term in Eq. (|26[) can be elaborated when higher moments 
are calculated. The latter are essential in order to find out a clear pattern. 

The figure Q] gives the results of the higher moments as calculated in HRG model. As introduced previously, 
raising the orders of moments results in new coefficients and new integrals. The earlier are partly characterized 
in Eq. (|26|) . The latter are estimated and depicted in the right panel of Fig. [TJ The integrals can be compared 
with the correlation functions (30[. We find that these integrals are related to (N) r , where r is the order of the 
moment. It is obvious that successive moments have a difference of about one order of magnitude. Therefore, 
higher order (N) r can be disregarded with reference to their vanishing contributions. The left panel depicts 
dimcnsionless non- normalized moments in dependence on the temperature T. 

In order to draw a picture about the contributions of the integrals appearing in Eqs. (|2"U)) - ([2~5j) . the phase 
space integral can be replaced by a series representation. For i-th particle, the pressure can be expressed as 

P(T, W ) = * T 4 £ (±)«+i e» * (n f ) n" 4 K 2 (n ^) , (27) 

n=l 

where the factor (±) n+1 represents fermions and bosons, respectively. In deriving previous expression, we 
expand the functions appearing in Eq. (|4|), for instance. The resulting expression is applicable as long as 
Hi < rrii is valid. Then, the derivatives with respect to fx/T leads to 

m r (T, W ) = * T 4 £ (±)" +1 e" * (» f ) ^ 4 K 2 (n f ) , (28) 

n=l 

which is valid for all moments of order r > 1. 
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For completeness, we mention here that taking into consideration one particle and its antiparticle is not rare 
in literature. In this case, the expression (|27p can be re-written as 

P(T, = * T 4 £ (±)" +1 (» Y ) K2 (" T) C ° Sh (" f ) ' (29) 

n— 1 

and accordingly, 

mrCT.A*,) = ^T 4 ^(ir 1 (n^) n^K 2 (n^) |Uosh („ |) . (30) 

n— 1 

Comparing Eq. (|28|) with Eq. ([27)) . we conclude that calculating higher moments is associated with a gradual 
change in n. Expression (|10p reflects almost the same behaviour, where higher moments are associated with 
change in the coefficients. 



B. Normalized higher moments 



The normalization of higher moments which can be deduced through derivatives with respect to the chemical 
potential [i of given charges, apparently gives additional insights about the properties of higher moments. From 
statistical point of view, the normalization is done with respect to the standard deviation a, which is be related 
to £. Therefore, it provides with a tool to relate moments with various orders to the experimental measurement. 
The susceptibility of the distribution give a measure for a. Should we are interested on the multiplicity, then 
the susceptibility is simply given by the derivative of first moments with respect to /i. It has been shown that 
the susceptibility is related to ~ £ 2 [if. The results of a in hadronic, bosonic and fermionic resonances are 
calculated at different fj, and given in Fig. [2] As discussed earlier, the strange quantum number has to remain 
conserved in high energy collisions. As per the standard model, this is one of the global symmetries in strong 
interactions. The procedure of keeping strange degrees of freedon conserved in HRG is introduced in Ref. 
This is the origin of the /i-dependence, especially in the bosonic gas. Although, baryo-chemical potential /i 
vanishes per definition, the chemical potential associated with strange quark /is remains finite. The dependence 
of /is on // is depicted in Fig. 1161 Another feature in these calculations is the assumption that the freeze-out 
boundary is determined by constant s/T 3 , where s is the entropy density [43| . At chemical freeze-out boundary, 
the dependence of a on /i is given in Fig. \5\ 
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Fig. 2: Results for a in hadronic (left), bosonic (middle) and fermionic (right) resonances are given in dependence on T 
for various baryo-chemical potentials (given in MeV). 

The "third" moment normalized to a 3 is known as 'skewness". For standard Gaussian distribution, the skew- 
ness is obviously vanishing. Therefore, the skewness is an ideal quantity probing the non-Gaussian fluctuation 
feature as expected near T c . The QCD CEP is conjectured to be sensitive to skewness. Experimentally, it has 
been shown that the skewness S is related to ~ £ 4 ' 5 [12]. The skewness for bosonic and fermionic resonance 
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£as, respectively, reads 



S b 



S, 



1 T 3 /2 / o °°C S ch[^] 4 Sinh[^i] P dk 
J o °°(Cosh[^] -l)" 1 k 2 dk 



[f °° (Cosh [^a] - 
, T 3/2 /o°° Csch[s^] 3 Sinh 
' ^ [fi> (Cosh [fl^] + 1 



: 3/2 



Sinh [ £ %^] 4 fc 2 rffc 



A: 2 dfc 



3/2 



(31) 



(32) 



At different chemical potentials, the skewness S is calculated in dependence on T and given in Fig. [3j It is 
worthwhile to notice that the values of S in the fermionic resonance gas is much larger than the S- values in the 
bosonic resonance gas, especially at small chemical potentials. But at large chemical potential, this observation 
is exactly opposite. The value of /z, at which the large (or small) skewness is flipped, seems to reflect the nature 
of phase transition, where the fermionic and/or bosonic degrees of freedom turns to be undistinguishablc. Such 
an observation shall be utilized in positioning the QCD CEP, section UTTl 




Fig. 3: At various baryo-chemical potentials (given in MeV), the skewness 5* for hadronic (left), bosonic (middle) and 
fermionic (right) resonance gas is given as function of T. 



The normalization of 4— th moment is known as heteroskedacity or kurtosis. It means varying volatility or 
more accurately, varying variance. Actually, the kurtosis is given by normalized 4— th moment minus 3. The 
subtraction of 3, which arises from the Gaussian distribution, is usually omitted [HI, 0, HEl . Therefore, the 
kurtosis is assumed to probe the non-Gaussian fluctuation feature as expected near T c and CEP. A sign change 
of skewness or kurtosis is conjectured to indicate that the system crosses the phase boundary [l6|, H(|. As HRG 
is valid below T c , the sign change is not accessible. It has been shown that kurtosis k is related to ~ £ 7 fl2l ). 



Kb 



T 

9i 



J °° {Cosh [a^et] + 2} Csch [^] 4 k 2 dk 



/ O °°(l-Coah[a^]) 1 k 2 dk 



3, 



tt 2 T 3 r °° {Cosh [a^] - 2} Sech fc 2 dk g 



/ °° (Cosh [STfii] + 1) 1 fc2 dk 



(33) 
(34) 



C. Products of higher moments 



There are several techniques to scale the correlation functions. The survey system's optional statistics module 
represents the most common technique i.e., Pearson or product- moment correlation. This module includes the 
so-called partial correlation which seems to be useful when the relationship between two variables is to be 
highlighted, while effect of one or two other variables can be removed. In present work, we study the products 
of higher moments of conserved quantities distributions. The justification of this step is that certain products 
can be directly connected to the corresponding susceptibilities in Lattice QCD simulation and related to long 
range correlations (l6l . ItTI . |48| . Seeking for simplicity, we start with the Boltzmann approximation. 
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1. Modified Bessel function: Boltzmann statistics 

When relativistic momentum integrals are replaced by summation over modified Bessel functions, the 
of higher moments in Boltzmann approximation read 

a 2 

(N) 



E oo =i(±r+ l e n^ 


(nf) 2 n 


~ 2 K 2 (n^) 


E - i(±) „+l e n^ 




-3K 2 (n^) ' 


E~i(±) n+1 e"^ 
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r 



k a 
~S~ 



2T 2 E^i(±)" +1 e"^ (n^) Z K 2 (n^) - 3 5 . 4 T 4 £~ 1 (±)«+ 1 e™^ (n^) n- 2 if 2 (n^ 



products 

(35) 
(36) 
), (37) 

(38) 



The origin of the second term in all expressions containing « is obvious. Expressions (|35[) - (|38[) justify the 
conclusions in [4^] that in Boltzmann approximation 

5(7 ~ 1, 



a 

(TV) 



while 



k a 



k a 

~~s 



l-3^exp 



(39) 
(40) 



These expressions are valid in the final state, which can be characterized by chemical and thermal freeze-out. 
In other words, these expressions are functions of the chemical potential. Using relativistic momentum integrals 
(see section Hi C 2\ shows that the products of moments result in a constant dependence on T [49j . 

2. Relativistic Momentum Integrals; quantum statistics 

The fluctuations of conserved quantities are assumed to be sensitive to the structure of the hadronic system 
in its final state. As mentioned above, crossing the phase boundary or passing through critical endpoint is 
associated with large fluctuations. Most proposed fluctuation of observables are variations of second order 
moments of the distribution, such as particle ratio [|| and charged dynamical measurement (50j . Then, the 
fluctuations are approximately related to £ 2 [l4| . 
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Fig. 4: The ratio a 2 /(N) is given in dependence on T at different chemical potentials fj, for hadronic (left), bosonic 
(middle) and fermionic (right) resonance gas. 



The ratio of standard deviation a 2 and the mean multiplicity (N) for fermions and bosons reads 

a 2 1 / °° (l±Csch[^]) _1 k 2 dk 



/"(lieV 4 ) k 2 dk 



(41) 
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where ± stands for fermions and bosons, respectively. The results are given in Fig. [4] We notice that the 
bosonic resonance gas results in smaller values than the fcrmionic one, especially in final state. Furthermore, 
we notice that a 2 / ' (N) decreases with increasing /i of the bosons at middle temperature. This is exactly the 
opposite in the fermionic resonance gas. This would give an explanation for the observation that the results 
in hadron resonances are not spread as in the other two sectors. Figure [5] shows the dependence of a 2 / (N) on 
H at the chemical freeze-out boundary, which is characterized by s/T 3 — 7. The ratio a 2 / (N) is equivalent to 
m^/m-i, Eq. (|3HI) - The dependence of bosonic and fermionic a on \i is given in Fig. [SJ as well. We notice that 
a is smaller than a 2 /(N), especially at small /i. At large \i both quantities are almost equal. 
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Fig. 5: The ratio a 2 /(N) and a are given in dependence on /i for boson and fermion resonances. 
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Fig. 6: The product Sa is given in dependence on T at different /if, values (given in MeV). 



The multiplication of skewness S by the standard deviation a is directly related to the thermodynamics of 
the number susceptibility of the lattice QCD. In HRG, the bosonic and fermionic products read 

1 /o -Csch[^] 4 Sinh[^] k 2 dk 

4 /^(l-Coshf™])" 1 k 2 dk ' 

, s f n °° Csch[^l 3 SinhT^l 4 k 2 dk t 
(Sa) f = 4^ L_iLJ L_^T_J _ 43 

f f o °°(Cosh[^] +1)" k 2 dk 

The product Sa is equivalent to TO3/7712, Eq. ()49j) . The results are given in Fig. [6j It is obvious that Sa ~ 1 
for either bosons or fermions. Then, for hadrons, Sa ~ 2. Nevertheless, the fine structure seems to reveal 
interesting features. 

The dependence of a and a 2 j (N) on T is illustrated in Figs. [2] and 01 respectively. It is obvious that 
both quantities have a monotonic behavior. Their dependences on (i are given in Fig. [5] Also, this type of 
dependences seems to be monotonic. As given in Fig. |6l the product S a has a characteristic dependence on 
T . Regardless, the tiny change, we notice that increasing fj, increases the bosonic S a product, but decreases 
the fermionic Sa product. In both cases, it forms a folding fan. The hadronic product makes an amazing 
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bundle in the middle (at a characteristic T). This will be discussed in details in section ITO1 When S is studied 
as a function of p, and given in Fig. 1151 we find that S(fi)b remains almost constant, while S(fi)f raises with 
increasing fi. We also notice that both curves cross at a certain point. This is not the case of the /i-dcpcndence 
of a and a 2 /(N), Fig. [5] So far, we conclude that starting from the third normalized moment, a non-monotonic 
behavior appears. Thus, the kurtosis and its products should be expected to highlight such a non-monotonic 
behavior 0, E3, [H| . 



The multiplication of kurtosis by a 2 called 52 1 is apparently equivalent to the ratio of 3-rd moment to 

2-nd order moment, Eq. ([50]) . In lattice QCD and QCD-like models, is found to diverge near the critical 
endpoint [IB, S3- I n HRG, the bosonic and fermionic products read 



2} Csch 



2T 



k 2 dk 



1 £ {Cosh gjg 

4 /^(i-Coshf^])" 1 k 2 dk 

1 J °° {Cosh - 2} Sech [^=^Y k 2 dk 

4~ / °° (Cosh [s^] + I)" 1 k 2 dk~ 



3jH_ J_ 
4tt 2 T 3 

3_£j 1 

4tt 2 T 3 



1 - Cosh 



T 



r H 






T 



A; 2 dfc, (44) 



k 2 dk. (45) 



When ignoring the constant term in Eqs. (f3"3"|) and (f3~4"| . then the second terms in the previous expressions 
entirely disappear. The results are given in Fig. [7] In the hadronic sector, the dependence of k a 2 on the 
temperature T at different /i-values is given in the left panel. We notice that increasing T is accompanied with 
a drastic declination in na 2 . Also, we hnd that no 2 flips its sign at large T. When comparing the thermal 
evolution of k a 2 with the one of S a, Fig. [51 we simply find that the latter is much drastically changing than 
the earlier. Also, when comparing their dependences on the chemical potentials at the freeze-out boundary, 
Figs. l9land[T4l it is apparent that the /^-dependence increases when the normalized fourth moment is included. 
The product kit 2 calculated at the freeze-out boundary leads to some interesting findings. First, no 2 almost 
vanishes or even flips its sign. Second, the T and /x corresponding to vanishing k a 2 are coincident with the 
phenomenologically measured freeze-out parameters. Third, the freeze-out boundaries of bosons and fermions 
are crossing at a point located very near to the one assumed by the lattice QCD calculations to be the QCD 
CEP, section [TV] 
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Fig. 7: The product k a 2 is given as a function of T at various fi- values for hadronic (left) , bosonic (middle) and fermionic 
resonance gas (right). 



The skewness and kurtosis can be combined through the standard deviation as follows. 

-\ir 2 T 3 r °° {2 + Cosh [^i] } Csch [™] 4 k 2 dk +\ 9i (j~ (l - Cosh k 2 dk 



k a 



\ it 2 T 3 r °° Csch [ £ %^f Sinh k2 dk 



/KG"* 



I ^ T 3 /o°° { Cosh - 2 ) Scch [^l 4 fc2 dk 



\9i 



/^(l + Coshf™])" 1 k 2 dk 



4 7T 2 T 3 J °° Csch [a^it] Sinh k 2 dk 



(46) 



(47) 



Also when ignoring the constant term in Eqs. Q33p and (|34]), the second terms in previous expressions disappear. 
The results of kct/S are given in Fig. [5] In the fermionic sector, na/S decreases with increasing T. Increasing 
H makes the decrease much faster. Although the drastic change in k<j/S, this behavior can be compared with 
Scr, Fig. |51 qualitatively. In the bosonic sector, there is a very slow decrease with increasing /i. The same 
sector in S a, Fig. |6l shows a very slow but an increasing dependence on T. Therefore, in the hadronic sector, 
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na/S is overall diminishing with increasing T. The declination becomes faster for larger \i. Fig. [9] illustrates 
the /^-dependence, i.e. na/S is estimated at the freeze-out curve. We notice that the bosons result in an almost 
unvarying na/S with growing fi. Therefore, the overall decrease in the hadronic sector is originated in the 
fcrmionic degrees of freedom. 




Fig. 8: ko/S as a function T at different chemical potentials for hadrons (left), bosons (middle) and fermions (right). 
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Fig. 9: The dependence of ko/S on the chemical potential /t. 
From Eqs. (|4Tj) -(l47 |) . we summarize that 



° - 7712 (AR) 

Sa^^, (49) 
m,2 

na 2 ee ^i, (50) 
na m 4 

"TT = — ■ (51) 
S m 3 



Obviously, other products would complete missing ratios, for example, 

Sx ^ ^, (52) 
mi 

m s (53) 

relating third and fourth moments to the first one, respectively. As discussed above, the susceptibility \ is 
equivalent to a 2 . These products seem to be sensitive to the volume independent multiplicities. 

The results of S x = S <r 2 are given in Fig. [TO] At the freeze-out boundary, the /i-dcpendcncc of S a 2 seems 
to be much stronger than the /(-dependence on S a, Fig. 1111 
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Fig. 11: The dependence of S a (left panel) and S a 2 (right panel) of bosons and fermions on the chemical potential fj, 
(given in MeV). 



The thermal evolution of k<j 4 /(N) is illustrated in Fig. Q21 Increasing T is accompanied by an almost overall 
increase of Ka i /(N) values. At low T, fluctuations appear. In bosonic sector, increasing [i leads to an obvious 
decrease in n a 4 /(N) so that the curves are distributed over a wide range. This is not the case for the fermionic 
sector. The left panel shows the results in the hadronic sector. At the freeze-out boundary, the values of 
kc a /(N) at various //-values are given in Fig. 1131 At large //, both bosonic and fermionic ratios vanish. With 
decreasing [i, there is a bounce. It reaching its maximum over ~ 0.25 GeV. Then both curves come closer to 
each others. They cross again over ~ 0.2 GeV. Continue decreasing // results in another bounce off. At smaller 
//, the fermionic sector switches again to positive values, while the bosonic sector results in large negative values 
of kct 4 /(A). For hadrons, na A /{N) remains almost negative. At very low /x, it turns to small positive values. 




Fig. 12: Thermal evolution of Ka i /(N) for hadrons (left), bosons (middle) and fermions (right) at different chemical 
potentials (given in MeV). 
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Fig. 13: The ratio na 4 /(N) are given in dependence on p for boson and fermion resonances. 



III. CHEMICAL FREEZE-OUT 



In rest frame of produced particle, the hadronic matter can be determined by constant degrees of freedom, 
for instance, s/T 3 (4/tt 2 ) = const [5j|. The chemical freeze-out is related to the particle creation. Therefore, 
the abundances of different particle species are controlled by chemical potential, which obviously depends on T. 
With the beam energy, T is increasing, while the baryon densities at mid-rapidity is decreasing. The estimation 
of the macroscopic parameters of the chemical freeze-out can be extracted from particle ratio. These parameters 
collected over the last two decades seem to fellow regular patterns as the beam energy increases [lo, HH . The 
higher moments have been suggested to control the chemical freeze-out, so that several conditions have been 
proposed [54| . 

As introduced in section Hi C 2i the thermal evolution of na 2 declines. This result seems to update previous 
ones [32l. [551] , where no 1 is assumed to remain finite and positive with increasing a. In present work, we find 
that the sign of no 1 is flipped at high T [56[. It has been found that the T- and /it-par ameters, at which the 
sign is flipped are amazingly coincide with the ones of the chemical freeze-out, Fig. Q31 Also, we find that 
the freeze-out boundaries of bosons and fermions are crossing at one point located at the hadronic curve. This 
point is close to the one that lattice QCD calculations suggest for QCD CEP, section ITVl The results are given 
in Fig. [14] Vanishing k <j 2 leads to 
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The rhs and lhs in both expressions can be re-written as 



1 - Cosh 



£i - Hi 



T 



( Cosh 
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T 



k 2 dk 



1 k 2 dk 



(54) 



(55) 



16— T 3 k 



9i 



-T 3 X 2 , 



which is valid for bosons and fermions. Then, the chemical freeze-out is fulfilled, if 

k{T,(i) = 3x 2 (7 1 , ) u). 



(56) 



(57) 



At chemical freeze-out curve, a naive estimation leads to £ ~ 3 1 / 3 fm. In doing this, it is assumed that the 
proportionality coefficients of k ~ £ 7 and x ~ £ 2 i are equal. In the heavy- ion collisions, £ has been measured [EH ]. 
Near a critical point, the experimental value ~ 2 — 3 fm (only factor 3 larger) agrees well with our estimation. 

At the chemical freeze-out curve, the intensive parameters T and (i which are related to the extensive proper- 
ties entropy and particle number, respectively, have to be determined over a wide range of beam energies. Left 
panel of Fig. [TJ] collects a large experimental data set. For a recent review, we refer to [2(| and the references 
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therein. The double-dotted curve represents a set of T and /j,, at which no 2 vanishes. It is obvious that this 
curve reproduces very well the experimental data. As given above, at this curve the normalized fourth moment 
K is equal to three times the squared susceptibility X- This new condition seems to guarantee the condition 
introduced in [43j : s/T 3 = const., right panel of Fig. [Mj at least over the range < u < 0.6 GcV. When 
excluding all degrees of freedom but the fermionic ones, s/T 3 decreases with increasing a (dashed curve). An 
opposite dependence is accompanied with the bosonic degrees of freedom (dotted curve). 
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Fig. 14: Left panel illustrates the chemical freeze-out boundary. The experimental data are given by the solid circles. 
The curves represent the results of hadron, boson and fermion gas determined when the sign of /tcr 2 is flipped. Right 
panel gives the corresponding s/T 3 values. The lattice QCD endpoint is marked by solid square. 



IV. QCD CRITICAL ENDPOINT 

In right panel of Fig. 1 141 it is interesting to notice that both fermionic and bosonic curves intersect at the 
hadronic curve at one point. Also, in the freeze-out diagram, all curves crossing at one point, left panel of Fig. 
[l~4l Furthermore, based on non-perturbative convergence radius, the critical endpoint as calculated in lattice 
QCD 0, Q given by solid square in Fig. Q3]is located very near to the crossing point. This is the motivation of 
Fig. M 
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Fig. 15: Left panel: Skewness and kurtosis at the freeze-out boundary for fermions and bosons. The curves of both 
quantities are crossing at almost the same fi-value. The corresponding T-values are almost equal. 

In Fig. 1151 the skewness and kurtosis of bosons and fermions calculated at the freeze-out curve are given in 
dependence on u (left panel) and center-of-mass energy yfs (right panel). It seems that the evolution of fermionic 
and bosonic skewness and kurtosis coincide at one point, marked with the vertical band. This fact would reflect 
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the nature of the phase transition. It would be the critical endpoint connecting cross-over with the first-order 
deconfinemcnt phase transitions. At QCD CEP, the phase transition is conjectured to be of second-order. It is 
worthwhile to mention the crossing point is amazingly coincident with the QCD CEP measured in lattice QCD 
[SHIi Fig. [Ml regardless its uncertainties. 

From Eqs. (|5T]) - (|32p and (|3"3"|) - (|3~i)) . following expressions have to be solved in /Lt, individually and/or depen- 
dency, in order to determine fi of the crossing point: 
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(58) 



(59) 



Due to the mathematical difficulties in dealing with these expressions, the integral over phase space has to be 
simplified. A suitable simplification is given in section III C II 



P(T, fig, Us) 
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1 ^ ~ (±)« 
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K 2 (n^ 
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(3n b + n s )fj, q - n s ^ s 
' T 



(60) 



where rif, and n a being baryon (strange) quantum number and \i q (/i s ) is the baryo- (strange) chemical potential 
of light and strange quarks, respectively The quarks chemistry is introduced in section [Ilj Accordingly, the 
difference between baryons and fcrmions is originated in the exponential function. For simplicity, we consider 
one fermion and one boson particle. Then, the baryo-chemical potential \i at the chemical freeze-out curve, at 
which the fermionic and bosonic skewness (or kurtosis) curves of these two particles cross with each other can 
be given as 



fib = iribHq = T In 



g b m 2 K 2 (^) 
g f m) K 2 (^f) 



(61) 



In the rclativistic limit, Ka(m/T) « 2T 2 /m 2 — 1/2 while in the non-relativistic limit Ka(m/T) « 
y/ 7rT/2mexp(— m/T)(l + 15T/8m). It is obvious that the bosonic and fermionic degrees of freedom play 
an essential role in determining Eq. (|61|) . Furthermore, it seems that the chemical potential of strange quark 
has no effect at the crossing point. 

The dependence of fxs on fj,b as calculated in HRG is given in Fig. 1161 As mentioned above, us is calculated 
to guarantee strange number conservation in heavy-ion collisions. At small /Lt&, fis has a linear dependence, 
Us = 0.25 Hb (Hooke's limit). At large fib, the dependence is no longer linear. 
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Fig. 16: us as a function of [it, in hadron resonance gas model (dotted curve). The linear fitting is given by solid curve. 



As discussed earlier, the critical behavior and the existence of QCD CEP can be identified by means of 
signatures sensitive to singular parts of the free energy, especially the ones reflecting dynamical fluctuations of 
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conserved charges, such as baryon number and charge density [581 ] . The reason that the fermionic and bosonic 
higher moments are crossing on a point very near the QCD CEP shall be elaborated in a forthcoming work [59j . 
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